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Abstract — Karst covers nearly half of Croatia's area. Croatian karst is part of the Dinaric 
karst and stretches from the Adriatic Sea to the Pannonian basin. Carbonate aquifers are the 
primary source of freshwater in karst part of Croatia and their protection is of the highest 
priority. Karst aquifers are sensitive to pollution since pollutants can easily enter the 
groundwater channel systems, often without prior self-purification process. Once pollutants 
enter the karstic underground they spread through the aquifer very quickly. Therefore, a 
thorough knowledge of the karst aquifer is essential for a timely and appropriate reaction to 
possible pollution incidents. Complexity of the karst landforms and groundwater networks 
requires implementation of a standard hydrogeological monitoring as well as unconventional 
methods of investigation. 

Analysis of stable water isotopes 7H and '8O proved to be helpful complementing method to 
a standard hydrogeological karst studies. We present the results of the stable isotope 
composition analysis of the coastal karst springs in the Bakar Bay and rain water collected in 
the hinterland of the springs. The local water supply company supervises the examined 
springs. During two years, spring water samples were collected on a weekly basis and rain 
samples were collected once a month. 

The stable isotope composition of the karst groundwater was modelled using autoregressive 
integrated moving average modelling. The study's main findings are: winter precipitation of 
Mediterranean origin dominates springs recharge, a dual porosity model that includes a 
fissure-porous aquifer and karstic channels is a fit for the studied systems, and autocorrelation 
function analysis revealed varying degrees of karstification in the hinterland of the studied 
springs. 


Introduction 


Karst aquifers are important sources of potable water not only in the Mediterranean 
region, but also globally [18]. They are recognized by their heterogeneous physical properties 
and complex flow patterns, which make investigation and description of their functioning 
challenging [5,14]. Many physicochemical parameters of groundwater could be used to 
characterize the hydrological behaviour of karst watersheds. However, in the case of non- 
stationary hydrological conditions, hydrochemical data should be supplemented with 
additional data such as environmental isotopes (e.g. water isotopes '8O and °H) to provide a 
reliable interpretation of karst basin functioning [15]. Some applications of water stable 


Referee List (DOI 10.36253/fup_referee_list) 
FUP Best Practice in Scholarly Publishing (DOI 10.36253/fup_best_practice) 
Diana Mance, Danijela Lenac, Maja Radišić, Davor Mance, Josip Rubinić, The use of 2H and 180 isotopes in 


the study of coastal karstic aquifers, pp. 525-534 © 2022 Author(s), CC BY-NC-SA 4.0, 10.36253/979-12-215- 
0030-1.48 


isotopes in karst aquifer research include the analysis of recharge processes [8,13] and water 
reservoir mixing [12], as well as the determination of residence times [19] and mean recharge 
elevations [16]. 

Although stable isotopes are typically used to supplement traditional hydrological 
methods [6], it has been demonstrated that when conventional parameters are unavailable, a 
thorough statistical analysis of oxygen and hydrogen stable isotope time series could be used 
for description of karst system hydrological behaviour [7]. Stable isotope content of Bakar 
Bay springs (Perilo - PER, Dobra - DB, and Dobrica - DBC) are discussed in the paper. Since 
there is no data on discharges for these springs, in our analysis we had to rely entirely on 
stable isotopes and statistical modelling. 


Materials and Methods 


An element's isotopes share the same number of protons but differ in the number of 
neutrons in their atomic nuclei. There are two stable isotopes of hydrogen: 'H and °H. The 
lighter hydrogen isotope accounts for ~ 99.985 % of total stable hydrogen, with the heavier 
isotope 7H accounting for the remaining ~0.015 %. There are three stable forms of oxygen: 
160, 70, and '80. The lightest one, is the most abundant, while !’O and '8O are less common 
in nature. 

There are nine different stable water configurations, the most common of which are 
'W'H'6O, 'H?H'°O, and 'H'H'8O. The masses of various stable water configurations differ, 
as do their physical and chemical properties. These differences result in isotopic fractionation 
or changes in stable isotope abundances at the beginning and end of physical, chemical, or 
biological processes. Because of fractionation, stable isotopes are sometimes referred to as 
"fingerprints" used in determining the origin of water [20]. 

The stable isotope composition of water is represented by 5!8O and 57H, with ô- 
value defined as the ratio of heavier to lighter isotope abundance in the sample (Rsampie) and 
the standard (Retandard): 5(%0) = Rsamte / Rstandard — 1. To express the stable isotope composition 
of water, the international VSMOW (Vienna Standard Mean Ocean Water) standard is used. 
Fresh water 5-values are typically negative, indicating a decrease in 7H or/and !8O abundance 
compared to the standard. 

Evaporation and condensation have a significant impact on the isotopic composition 
of water (Fig. 1). Water vapour is depleted in heavy isotopes when compared to the evaporating 
body. In contrast, rain contains more heavy isotopes than residual vapour. Air temperature 
has a strong influence on precipitation isotope content, resulting in higher -values in the 
summer and lower 5-values in the winter [11]. There is a linear correlation between °H and 
180 values in natural waters that are not affected by evaporation (67H = 8-6'8O + 10 %o), 
known as a Global Meteoric Water Line (GMWL) [2]. The calculation of regression lines for 
local precipitation (Local Meteoric Water Line - LMWL) and ground water (Local 
Groundwater Line - LGWL) is a standard procedure in isotope hydrology. If there is no 
further evaporation, the isotopic composition of the precipitation usually remains unchanged 
once it enters the underground [11]. Nonetheless, mixing of different water masses causes 
changes in the stable isotope composition of groundwater [10]. 
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Figure 1 — The atmospheric part of water cycle and changes in stable isotope content 
caused by water evaporation and condensation. 


In 1964, Dansgaard introduced the concept of d-excess, which is defined as follows: 
d-excess = 57H - 8-5!8O [3]. D-excess values of around 10 %o indicate precipitation originating 
in the Atlantic Ocean, whereas values greater than 15 %o indicate precipitation originating in 
the Mediterranean [4]. 

The study area is the Bakar Bay spring discharge zone in western Croatia (Fig. 2). 

For two years, weekly spring water samples were collected. We used a set of 
precipitation isotope data from the Kukuljanovo (KUK), Skalnica (SKAL), and Platak 
(PLAT) rain gauging stations to compare the isotopic content of precipitation and 
groundwater (Fig. 2). Some of these data have already been used in hydrogeological analyses 
of the Rjecina River chacment [7, 8]. To prevent evaporation, monthly totals were collected 
in 3.5-litre rain gauges containing 100 ml of paraffin oil. After being separated from the oil, 
the precipitation samples, same as groundwater samples, were stored in double-capped 
HDPE bottles. A Delta"SXP (Thermo Finnigan) isotope ratio mass spectrometer, with an 
HDOeq48/24 (IsoCal) equilibration unit and a Dual Inlet system (Thermo Finnigan) as 
peripherals, was used for the stable isotope measurements. For 5!8O, measurement precision 
was better than 0.1 %o, and for 87H, it was better than 1 %o. 
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Figure 2 — The main map: a hydrogeological map of the study area showing basins of 
Rjecina River and coastal springs in the Bakar Bay. The map also shows: the main 
karst springs (RJ — Rjecina Spring, ZV — Zvir, MWs — Martinšćica wells, PER — 
Perilo, DB — Dobra, DBC — Dobrica), the Rjecina River (blue line on the main map), 
and rain gauge stations for sampling cumulative monthly precipitation (KUK — 
Kukuljanovo, SKAL — Škalnica, PLAT — Platak). Upper right corner: the position of 
the study area in Europe. 
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To analyse the groundwater stable isotope series, we used univariate time series 
analysis methods such as autocorrelation function (ACF) and autoregressive integrated 
moving average modelling (ARIMA) [1]. We used graphical representations of ACF to draw 
a conclusions about aquifer's karstification and retention capabilities, as well as to select the 
appropriate ARIMA model. Prior to ARIMA, we used an additive component model (6 =/+ 
s t+sc) to test the groundwater stable isotope time series for trends and seasonal variations, 
where / is a linear trend component determined by linear regression, s is a seasonal 
component determined by periodic regression, and sc is a stochastic component. The results 
were interpreted using the 0.05 level of statistical significance. 


Results and Discussion 


For the study area and period, LMWLs for the warm part of the hydrological year 
(8°H = 8.05:8'8O + 8.52 %o; R°=0.89) and cold part of the hydrological year (57H = 8.04-5!8O 
+ 13.96 %o; R?=0.99) were already presented in [7]. The division of the LMWL into summer 
and winter seasons was justified by previous studies. These have shown that the rainy season 
from October to April contributes more to groundwater recharge in the study area than the 
precipitation that falls during the summer season [7,8]. Figure 3 shows that the -values of 
all three springs (PER, DB, and DBC) agree well with the LMWL for the cold season, 
indicating that the groundwater is primarily fed by winter precipitation. LGWL for PER, DB, 
and DBC, as shown in Fig. 3 (67H = 8.37:'8O + 17.6 %o; R=0.95), further supports this 
conclusion. 

Figure 4 shows large variations in the d-excess for rainwater: values in the winter 
months are significantly higher than values in the summer months of the hydrologic year. 
Since precipitation is generally lower in the summer months than in the winter months, a 
weighted average of d-excess was calculated: KUK (12.9 %o), SKAL (14.71 %o) and PLAT 
(14.66 %0) [7]. Values for SKAL and PLAT correspond to the average values of groundwater 
d-excess for DBC (14.67 + 0.73) %o, DB (14.82 + 0.73) %o, and PER (14.53 + 0.73) %o. 

Due to the strong and statistically significant correlation between 8!°O and 8H 
values (R=0.97, p < 0.001), only 5'8O values were used for time series analysis. The ACF of 
the groundwater 5'8O series was analyzed to obtain information about the system's 
karstification. Figure 5 shows the ACFs of the 5!8O series for PER, DB, and DBC. 

The shape of the ACF discharge graph is commonly used to determine the retention 
capacity and degree of karstification of the karst aquifer [17]. The karst system's “memory 
effect” is the time required for auto-correlation coefficients r(k)xx to fall below 0.2 [9]. We 
used the same reasoning to interpret 5'3O ACFs (Fig. 5). The ACFs do not decrease 
continuously but exhibit two types of decline, supporting the interpretation of the aquifers’ 
dual nature. Initially, all three ACFs drop rapidly with steep slopes. The gradients begin to 
decline after five weeks. The ACF for PER decreases more rapidly, and its autocorrelation 
factors become statistically insignificant after eight weeks; whereas the ACFs for DB and 
DBC decrease more gradually, and their autocorrelation factors remain statistically 
significant for eleven weeks. The slower declines in DB and DBC ACFs suggest that DB and 
DBC have a greater ability to retain water than PER. This could also indicate a more 
developed karst channel system in the case of PER versus DB and DBC. 
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Figure 3 — Correlation diagram of ôH and 8!8O values for precipitation collected at 
Kukuljanovo (KUK), Skalnica (SKAL), and Platak (PLAT) stations, as well as groundwater 
collected at Dobrica (DBC), Dobra (DB) and Perilo (PER). Local meteoric water lines 
(LMWLs) for precipitation in the warm and cold seasons of the hydrological year, as well as 
a global meteoric water line (GMWL) are also presented. 
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Figure 4— D-excess values for precipitation collected at the stations Kukuljanovo (KUK), 
Skalnica (SKAL), and Platak (PLAT), and spring water from Dobrica (DBC), Dobra (DB) 
and Perilo (PER). 
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Figure 6 — Actual õ!ŝO data from Dobrica, fitted model (created with linear regression, 
periodic regression, and ARIMA (1,0,0)), and residuals. The circles in the upper part 
of the graph represent the measured data for 5!8O in water. A red line represents the 
fitted model. Deviations from the fitted model are represented by a black solid line in 


the lower part of the graph. 
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Since ARIMA results for PER and DB have already been published [8], this paper 
includes only DBC results. Significant linear and seasonal trends were observed: 5!8O = 
— 8.07 + 0.004-¢ — 0.32-cos((2m/52)-t+2.32)) (R?=0.65), where ¢ is a time period (a week). 
There were no seasonal changes in the 5!8O values of the monthly precipitation samples 
observed during the study period [7,8]. As a result, we conclude that seasonality in 
groundwater 5'8O time series is caused by aquifer characteristics such as rapid infiltration 
and discharge of heavy rain, as well as the predominance of winter precipitation in 
replenishing groundwater reserves. Indeed, the strongest 5!8O shifts to higher values were 
observed during the cold part of the hydrological year at the time of the heaviest precipitation. 
This indicates that precipitation infiltrated rapidly into the subsurface and reached springs in 
a short time. In contrast, groundwater had the lowest delta values during the warm part of the 
hydrologic year and these were without significant fluctuations. There were no heavy rain 
events during this period and we can conclude that the water sampled at the springs during 
this period came from the aquifer water reserves. The low delta values indicate that snowmelt 
probably plays the most important role in recharging the aquifer. This situation suggests that 
the systems analyzed are best described by the so-called dual porosity model. Such a model 
includes a fissure porous part of the aquifer characterized by baseflow whose activity is most 
evident during dry seasons. The second part of the model consists of a developed network of 
karst channels that respond rapidly to intense precipitation events. 

ARIMA modeling of the 8!8O series was performed after seasonal and linear trends 
were removed from the series and stationarity of the residual series was achieved. AR (1) 
(i.e., ARIMA (1,0,0) model) proved to be the best model for all three springs. The actual 
data, fitted data, and 5'8O DBC residuals are shown in Fig. 6. The largest deviation from the 
statistical models was found for samples collected during heavy rainfall. For PER they were 
up to 1.2 %o, while for DB and DBC they were not higher than 0.5 %o. This difference in 
residuals could be explained by the greater karstification of the hinterland of PER compared 
to DB and DBC. The good agreement of measured and modelled values (R? = 0.81) obtained 
by ARIMA (Fig. 6), also supports the idea of the DBC hinterland having a lower 
karstification degree and higher retention capability than PER (R? = 0.76). 


Conclusions 


To protect drinking water sources and use them sustainably, we need to understand 
how aquifers work. Karst aquifers are well-known for their complexity and research 
difficulty. Scientists have long agreed that multidisciplinary research is essential, and isotopic 
hydrology is an excellent example of multidisciplinary karst aquifer research. 

This paper presents the findings of a two-year sampling at three karst springs in 
Bakar Bay (Croatia), as well as rain gauge stations in their vicinity. Based on the isotopic 
composition of the collected samples, we concluded that these karst springs are primarily fed 
by winter precipitation. Temporal changes in groundwater isotopic composition show a 
seasonal pattern not seen in precipitation. As a consequence, the seasonal oscillation of 
isotopic composition is interpreted as a function of the karst system itself, i.e. as a dual 
porosity model consisting of a fissure-porous aquifer (characterized by baseflow during the 
dry season) and highly developed karstic channels (characterized by rapid infiltration and 
strong discharge during heavy rainfall). The analysis of auto-correlation functions and time 
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series show that there is a difference in the degree of karstification of individual springs, with 
a higher degree of karstification indicating a greater sensitivity to potential pollution. 
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